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ABSTRACT 

We report results from radiation hydrodynamical simulations of the collapse of molecular 
cloud cores to form protostars. The calculations follow the formation and evolution of the 
first hydrostatic core/disc, the collapse to form a stellar core, and effect of stellar core for- 
mation on the surrounding disc and envelope. Past barotropic calculations have shown that 
rapidly -rotating first cores evolve into 'pre-stellar discs' with radii up to ~ 100 AU that may 
last thousands of years before a stellar core forms. We investigate how the inclusion of a re- 
alistic equation of state and radiative transfer alters this behaviour, finding that the qualitative 
behaviour is similar, but that the pre-stellar discs may last 1.5 — 3 times longer in the more 
realistic calculations. The masses, radii, and lifetimes of the discs increase for initial molec- 
ular cloud cores with faster rotation rates. In the most extreme case we model, a pre-stellar 
disc with a mass of 0.22 and a radius of ^ 100 AU can form in a l-M© cloud and last 
several thousand years before a stellar core is formed. Such large, massive objects may be 
imaged using ALMA. Fragmentation of these massive discs may also provide an effective 
route to binary and multiple star formation, before radiative feedback from accretion onto the 
stellar core can inhibit fragmentation. Once collapse to form a stellar core occurs within the 
pre-stellar disc, the radiation hydrodynamical simulations produce qualitatively different be- 
haviour from the barotropic calculations due to the accretion energy released. This drives a 
shock wave through the circumstellar disc and launches a bipolar outflow even in the absence 
of magnetic fields. 

Key words: accretion, accretion discs - hydrodynamics - radiative transfer - stars: formation 
- stars: low-mass, brown dwarfs - stars: winds, outflows. 



1 INTRODUCTION 

The first numerical calculations of the collapse of a molecular cloud 
core to stellar core formation, and beyond, were performed by |Lar-| 
|son| ( |1969] ) more than four decades ago. These one-dimensional ra- 
diation hydrodynamical calculations revealed the main stages of 
protostar formation: an almost isothermal collapse until the inner 
regions become optically thick, the almost adiabatic formation of 
the first hydrostatic core (typical radius ^ 5 AU and initial mass 
^ 5 Mj), the growth of this core as it accreted from the infalling 
envelope, the second collapse within this core triggered by the dis- 
sociation of molecular hydrogen, the formation of the stellar core 
(initial radius ^ 2 R© and mass ^1.5 Mj), and, lastly, the long 
accretion phase of the stellar core to its final mass. More recent one- 
dimensional studies (e.g. |Masunaga & Inutsuka|2000[|Commergon| 
[et al.|201 1 ) have not altered this picture. 

However, in one-dimension, the effects of rotation and mag- 
netic fields can not be examined. |Larson| ( |1972| ) began per- 
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forming two-dimensional calculation of rotating clouds soon af- 
ter his ground-breaking one-dimensional calculations. But multi- 
dimensional calculations are considerably more time consum- 
ing and it was not until almost two decades after the first one- 
dimensional models that the first two dimensional calculations 
were able to follow the collapse to the formation of the stellar core 
( |Tschamuter 1 1 9 8 7|) . These calculations showed that with rotation, 
both the first hydrostatic core and the later-formed stellar core be- 
came rotationally-flattened and allowed the shock structures in the 
accretion flows to be studied in detail (e.g. |Tscharnuter et al.|2009| l. 

It was not until nearly three decades after |Larson| s original 
calculations that the first three-dimensional calculations following 
the collapse to stellar core formation were possible. |Bate| ( |1998| l 
performed three-dimensional calculations of rotating molecular 
cloud cores using a barotropic equation of state to mimic the effects 
of a realistic equation of state and radiative transfer. He found that 
if the original molecular cloud core was rotating rapidly enough, 
the rotationally-flattened first hydrostatic core could be dynami- 
cally unstable to the growth of non-axisymmetric perturbations. It 
deformed into a bar, followed by the wrapping up of the ends of the 
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bar into a disc with strong trailing spiral arms. Gravitational torques 
removed angular momentum and rotational support from the in- 
ner regions of the first core, quickening the onset of the second 
collapse and inhibiting fragmentation during the second collapse 
phase. The development of a bar-mode and spiral structure is ex- 
pected for rapidly-rotating polytropic-like structures (e.g.'Durisen" 
|et al.|1986t . Such instabilities occur when the ratio of the rotational 
energy to the magnitude of the gravitational potential energy of the 
first core exceeds /3 = 0.274. Bate was also the first to point out 
that because a rapidly-rotating first core develops into a disc be- 
fore the stellar core forms, the disc forms before the star. Rather 
than hydrostatic cores, such structures are better described as 'pre- 
stellar discs'. Once the second collapse occurred and a stellar core 
formed, ^Bate found an inner disc formed around the stellar core 
and grew in radius as material with larger angular momentum fell 
in. Bate speculated that if the calculation were able to be followed 
further, the outer radius of the inner disc would eventually grow to 
join on to the outer disc, leaving the stellar core surrounded by a 
single large disc. 

The next major advance in this field was to begin including the 
effects of magnetic fields. |Tbmisaka| ( [2002| ) was the first to follow 
the collapse of a molecular cloud core to stellar densities including 
magnetic fields. Using two-dimensional calculations, he found that 
both the first core and stellar cores launched magnetically-driven 
outflows. This work has been followed by a number of papers us- 
ing three-dimensional calculations to investigate magnetised col- 
lapse and outflows (Machida et al. 2005 2006 2008 J and the rota- 
tion rates and magnetic field strengths in the first and stellar cores 
( [Machida et al.|200"7] l. At the same time, the development of non- 
axisymmetric structure at the first core stage and the formation of 
a pre-stellar disc discovered by Bate (1998 ) and fragmentation has 
been studied further both withou t (^Saigo & Tomisaka^2006| |Saigo| 
|et al. |2008[ Saigo & Tomisaka| 201l| and with magne tic fields 
dMachida et al.|2005||2010[|Macliida & Matsumoto|201 With the 
increasing capability of observational instruments and, in particu- 
lar, looking forward to the capabilities of the Atacama Large Mil- 
limeter/submillimeter Array (ALMA), some of these works also 
began to investigate the observational signatures of the first core 
phase such as their luminosities ( Saigo & Tomisaka'2006 ), spectral 
energy distributions (SEDs), and images ( Saigo & Tomisaka 201 11. 

However, all three-dimensional studied listed above have used 
barotropic equations of state rather than realistic equations of state 
with radiative transfer. The first three-dimensional calculations in- 
cluding radiative transfer that followed collapse to the point of 



stellar core formation (but not beyond) were W hitehouse & Bate 
(|2006|), using t he flux-limited diffusion approximation, and |Sta- 
jmatellos et al.| ( |2007| ), using a radiative cooling approximation. 
^Bate (201 Of investigated the evolution of both first and stellar cores 
using radiation hydrodynamical calculations. He found that the 
evolution up until stellar core formation was qualitatively similar 
to that obtained with a barotropic equation of state, including bar- 
mode instabilities and stellar core formation. However, following 
the formation of the stellar core, the energy released has a dramatic 
effect on the surrounding disc and envelope and launches a tempo- 
rary outflow even in the absence of a magnetic field. Recent two- 
dimensional radiation hydrodynamical simulations that follow the 
collapse to stellar densities and well beyond as the stellar core ac- 
cretes from the inf ailing envelope (Schonke & Tscharnuter|20l"T] ) 
show a similar effect, where the formation of the stellar core drives 
a temporary outflow. 

The most complete three-dimensional calculations of the col- 
lapse of molecular cloud cores to date include radiation hydrody- 



namics and magnetic fields ( Tomida et al.||201()| i; [Tomida et al.| 
|2010] 3; [Commergon et al.||2011| ). As with the past magnetised 
barotropic calculations, the first core is found to launch outflows, 
but now the thermodynamics of the gas is more properly treated. 
However, these models only study the evolution of the first core, 
and do not follow the second collapse and formation of the stellar 



In this paper, we extend the study of Bate (2010) who pre- 
sented the first three-dimensional radiation hydrodynamical calcu- 
lations to follow the collapse beyond the formation of the stellar 
core. Bate only presented results for one particular value of the ro- 
tation rate of the molecular cloud core. Here we investigate how the 
evolution of the first core or pre-stellar disc varies with different 
rotation rates, and we compare the results obtained using radiation 
hydrodynamics with those obtained using a barotropic equation of 
state. We also study how the thermally-launched outflows follow- 
ing stellar core formation depend on the rotation rate, and how the 
discs evolve. Finally, we perform many of the calculations at sev- 
eral different numerical resolutions to test the effects of numerical 
resolution and examine convergence. 



2 COMPUTATIONAL METHOD 

The calculations presented here were performed using a three- 
dimensional smoothed particle hydrodynamics (SPH) code based 
on the original version of Benz (1990, Benz et al.| 19901), but sub- 
stantially modified as described in Bate et al. ( 1995 ), Whiteh ouse, 
Bate & Monag hanI ( [20051 ), [Whitehouse & Bate (2006 ), Pric~ 
Bate| ([2007 ), and parallelised using both OpenMP and MPL 

Gravitational forces between particles and a particle's near- 
est neighbours are calculated using a binary tree. The smoothing 
lengths of particles are variable in time and space, set iteratively 
such that the smoothing length of each particle h — 1.2{m/ pY^^ 
where m and p are the SPH particle's mass and density, respectively 
(see Price & Monaghan|2007] for further details). The SPH equa- 
tions are integrated using a second-order Runge-Kutta-Fehlberg in- 
tegrator with individual time steps for each particle ( [Bate et al 
|1995|). To r educe numerical shear viscosity, we use the [Morris & 
|Monaghan[ \ 1 997| ) artificial viscosity with varying between 0.1 
and 1 while /3v = (see also |Price & Monaghan|2005| ). 



2.1 Equation of state and radiative transfer 

The calculations presented in this paper, use two different types of 
equation of state. The first is a barotropic equation of state, almost 
identical to that used by Bate ( 1998 ), where the temperature of the 
gas depends only on its density. In this case, the thermal pressure 
of the gas p = Kp^ . The value of the effective polytropic exponent 
?7, varies with density as 



1, 1.0x10"^^ 

7/5, 1.0 X 10"^^ < p ^ 5.7 X 10"^ 

1.15, 5.7 X 10"^ < />< 1.0 X 10"^ 

5/3, />>1.0xl0"^ 



(1) 



We take the mean molecular weight of the gas to be /i = 2.38. The 
value of K is defined such that when the gas is isothermal K — cl, 
with the sound speed Cg = 2.04 x 10^ cm s~^, and the pressure is 
continuous when the value of rj changes. 

The second type of calculation is performed using radiation 
hydrodynamics. In this case, we use an ideal gas equation of state 
p = pTlZ/ p, where T is the gas temperature, and IZ is the gas 
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Figure 1. The time evolution of the maximum density (left panels) and gas temperature (right panels) during the barotropic calculations (top panels) and the 
radiation hydrodynamical calculations (bottom panels) of the collapse of molecular cloud cores with different initial rotation rates. For each panel, the different 
lines are for cloud cores with /3 = 0, 5 x IQ-^, 0.001, 0.005, 0.01, 0.04 from left to right. The free-fall time of the initial cloud core, % = 1.8 x 10^^ ^ 
(56,500 yrs). Each calculation was performed with 10^ SPH particles, except the barotropic calculation with (3 = 0.04 which used 3 x 10^ particles. 



constant. The thermal evolution takes into account the translational, 
rotational, and vibrational degrees of freedom of molecular hydro- 
gen (assuming a 3:1 mix of ortho- and para-hydrogen; see |Boley] 
|et al.|2007] ). We also include molecular hydrogen dissociation, and 
the ionisations of hydrogen and helium. The hydrogen and helium 
mass fractions are X = 0.70 and Y = 0.28, respectively. The 
contribution of metals to the equation of state and the thermal evo- 
lution is neglected. Two temperature (gas and radiation) radiative 
transfer in the flux-limited diffusion approximation is implemented 
as described by Whitehouse et al.| ( |200 5 ) and Whitehouse & Bate 
( [200 6 ), except that the standard explicit SPH contributions to the 
gas energy equation due to the work and artificial viscosity are 
used when solving the (semi-)implicit energy equations to provide 
better energy conservation. We assume solar metallicity gas, using 
the interstellar grain opacity tables of Pollack et al.| ( |1985) and the 
gas opacity tables of Alexander, ( ,1975i ) (the IVa King model) (see 
[Whitehouse & Bate|2006l L 



2.2 Initial conditions 

The initial conditions for the calculations are identical to those of 
|Batel ( [T998] ) and [Bate] pOTO] ). We follow the collapse of initially 
uniform-density, uniform-rotating, molecular cloud cores of mass 
M = 1 M© and radius R = 7 x 10^^ cm. The free-fall time 
of the initial clouds is ts = 1.785 x 10^^ s (56,500 yrs). The 
ratios of the thermal and rotational energies to the magnitude of 
the gravitational potential energy are a = 0.54 and /3, respectively. 
We have performed calculations with /3 = (not rotating), and /3 = 
5 X 10~^, 0.001, 0.005, 0.01 and 0.04. Typically, each calculation 
was performed twice, once using the barotropic equation of state, 
and again using radiation hydrodynamics. 

To satisfy the resolution criterion of [Bate & Burkert| ( [1997| ) 
that the minimum Jeans mass during the calculation contains at 
least ^ 2A^neigh = 100 particles, we require at least 1x10^ equal- 
mass particles. To test for convergence, we performed calculations 
using 1 X 10^ 3 X 10^ 1 x 10^ and 3 x 10^ equal-mass SPH 
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Table 1. A summary of the type of calculations performed, and their nu- 
merical resolutions. From left to right, the columns give the type of initial 
condition (i.e. the values of /3, the ratio of the rotational energy to the mag- 
nitude of the gravitational potential energy for the molecular cloud cores), 
the equations of state used (i.e. barotropic equation of state, or radiation hy- 
drodynamics with a realistic equation of state), and which numbers of SPH 
particles were used to perform calculations. 

particles. Only the /3 = 0.005 radiation hydrodynamical calcula- 
tion was performed with all these resolutions, but most calculations 
were performed with at least two different resolutions (see Table[TJ. 
Calculations using the highest resolution of 3 x 10^ SPH particles 
were only performed for the /3 = 0.001 cases and radiation hydro- 
dynamical 13 — 0.005 case. In this paper, unless otherwise stated, 
when a particular set of initial conditions and equation of state is 
discussed or a figure is presented, the calculation performed with 
the highest resolution is used. We discuss the degree to which the 
simulations are numerically converged in Appendix A and at vari- 
ous points in Section [3] 

The calculations were performed on the University of Exeter 
Supercomputer, an SGI Altix ICE 8200. 



3 RESULTS 

The collapse of each molecular cloud core up until the formation of 
the stellar core proceeds in a manner that is qualitatively similar to 
those reported from previous three-dimensional calculations with- 
out magnetic fields and using barotropic equations o f state (jBate 
|1998,,Saigo & Tomisaka|[2006, .Saigo et al.^2008, Machida et al. 



This is mimicked in the barotropic calculations by the change in 
the value of 77 from 1 to 7/5 (equation[TJ. This increasing tempera- 
ture leads to the formation of a pressure-supported 'first hydrostatic 
core' ( Larson 1 969 ), which can be seen in Fig. [T] when the initial 
collapse stalls with central (maximum) densities ^ 10~^^ g cm~^ 
and temperatures of ^ 30 — 120 K, depending on the degree of 
rotational support (i.e. cores that rotate more quickly have lower 
maximum temperatures). Again, apart from the small offset in the 
time of first core formation, the barotropic and radiation hydrody- 
namical calculations give very similar results to this point. 

Without rotation (/3 = 0), the first core has an initial mass 
of 5 Jupiter masses (Mj) and a radius of 5 AU (in agreement 
with |Larson 1969 ). However, with higher initial rotation rates of the 
molecular cloud core, the first cores become progressively more 
oblate (Figures [2] |3] and|4]). For example, with /3 = 0.005 using 
radiation hydrodynamics, before the onset of dynamical instability, 
the first core has a radius of 20 AU and a major to minor axis 
ratio of «4:1 (first panel in each third row of Figs. |3] and |4]). With 
13 — 0.01, the first core has a radius of ^ 30 AU and a major to 
minor axis ratio of ~6:1 (fourth rows of these figures). Thus, for 
the higher rotation rates, the first core is actually a pre-stellar disc, 
without a central object. As pointed out by |Bate| (1998| ), [Machidaj 
|et al.|p010| ), and |Bate| ( [20T0l ), the disc actually forms before the 
star. For the very highest rotation rates {(3 — 0.04), the first core 
actually takes the form of a torus or ring (first panel in each bottom 
row of Figs. |3] and |4| in which the central density is lower than the 
maximum density. 

3.1 Slowly-rotating first hydrostatic cores 

The evolution of the first core up until the point of stellar core for- 
mation depends on its rotation rate. Non-rotating and slowly rotat- 
ing cores evolve as they accrete mass from the surrounding infalling 
envelope with their central densities and temperatures increasing 
(Figs. [1] and |5] calculations with (3 ^ 0.001). In the radiation hy- 
drodynamical calculations, when the central temperature exceeds 
^ 2000 K, molecular hydrogen begins to dissociate, leading to a 
second hydrodynamic collapse deep within the first core ( |Larson| 
|I969| ). The formation of the stellar core occurs just a few years af- 
ter the onset of the second collapse, during which the maximum 
density increases from - 10"^ to ^ 0.1 2 cm and the max- 
imum temperature increases from ^ 2000 to > 60, 000 K. The 
stellar core is formed with a mass of Mgc ~ 1.5 Mj and a radius of 
Rsc ~ 2 R© . Without rotation, the stellar core accretes the remnant 
of the first core in which it is embedded in 10 years and then ac- 
cretes the envelope ( |Larson|1969| ), though with three-dimensional 
calculations we only follow the calculations for 50 — 100 years 
after stellar core formation. 

Although these stages are qualitatively the same in the 
barotropic calcuations, there is a much greater difference in the 
evolution of the first core between the barotropic and radiation hy- 
drodynamical calculations than for the phase of the collapse prior 
to first core formation. To make this clear, in Fig. |5] the time evo- 
lution of maximum density and maximum temperature during the 
calculations is replotted with t = set to the time of stellar core 
formation (defined as the time when the maximum density reaches 
10~^ g cm~^). This allows us to clearly see the amount of time 
spent between first core formation and stellar core formation in 
each of the calculations. The barotropic results are plotted using 
dashed lines, while the solid lines give the radiation hydrodynam- 
ical results. The evolution time of the first core or pre-stellar disc 
is longer using radiation hydrodynamics than using the barotropic 



2010[|Machida & Matsumoto|2011| ) for both the barotropic and ra- 
diation hydrodynamical calculations presented here. 

Fig. [1] gives the evolution with time of the maximum den- 
sity and temperature for clouds with different initial rotation rates, 
with each calculation performed using 10^ particles (except for the 
baratropic /3 = 0.04 case). The initial collapse is isothermal until 
the maximum density exceeds 10~^^ g cm~^ using the barotropic 
equation of state. Using radiation hydrodynamics the evolutions 
to this density are also almost isothermal. However, slight heat- 
ing when using radiation hydrodynamics does slow the collapse 
marginally, leading to the times taken to form the first hydrostatic 
cores being slightly longer 0.01 — 0.015 tff) than in the corre- 
sponding barotropic calculations. 

In the radiation hydrodynamical calculations, as the heating 
rate in the central regions exceeds the rate at which the gas can 
cool, the gas begins to heat up and the collapse enters an almost adi- 
abatic phase where the temperature rises as the gas is compressed. 
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Figure 2. Snapshots of the column density viewed parallel to the rotation axis during the evolution of the barotropic calculations of the collapse of molecular 
cloud cores with different initial rotation rates. From top to bottom, the different rows are for cloud cores with /3 = 5 x 10~^, 0.001, 0.005, 0.01, 0.04. Note 
that the spatial scale is different for each row, with each panel measuring 1600 y/^^ AU across (i.e. from 36 to 320 AU). The free-fall time of the initial cloud 
core, tff = 1.8 x 10^^ s (56,500 yrs). Each calculation was performed with 10^ SPH particles, except for the /S = 0.001 case which used 3 x 10^ SPH 
particles and the /5 = 0.04 case which used 3 x 10^ particles. The evolution is almost identical when using 3 x 10^ particles or more, but the latter are 
slightly more detailed. 
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Figure 3. Snapshots of the column density viewed parallel to the rotation axis during the evolution of the radiation hydrodynamical calculations of 
the collapse of molecular cloud cores with different initial rotation rates. From top to bottom, the different rows are for cloud cores with /3 = 5 x 
10~^, 0.001, 0.005, 0.01, 0.04. Note that the spatial scale is different for each row, with each panel measuring 1600 AU across (i.e. from 36 to 320 
AU). The free-fall time of the initial cloud core, tff = 1.8 x 10^^ s (56,500 yrs). Each calculation was performed with 10^ SPH particles, except for the 
(3 = 0.001 and /3 = 0.005 cases which used 3 x 10^ SPH particles. The evolution is almost identical when using 10^ or 3 x 10^ particles, but the latter are 
slightly more detailed (see Appendix A). 
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Figure 4. Snapshots of the column density viewed perpendicular to the rotation axis during the evolution of the radiation hydrodynamical calculations 
of the collapse of molecular cloud cores with different initial rotation rates. From top to bottom, the different rows are for cloud cores with /3 = 5 x 
IQ-^, 0.001, 0.005, 0.01, 0.04. Note that the spatial scale is different for each row, with each panel measuring 1600 AU across (i.e. 36, 50, 1 14, 160, or 
320 AU). The free-fall time of the initial cloud core, fff = 1.8 x 10^^ s (56,500 yrs). Each calculation was performed with 10^ SPH particles, except for the 
(3 = 0.001 and /3 = 0.005 cases which used 3 x 10^ SPH particles. The evolution is almost identical when using 10^ or 3 x 10^ particles, but the latter are 
slightly more detailed (see Appendix A). 
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Figure 5. The time evolution of the maximum density (upper panels) and gas temperature (lower panels) during the radiation hydrodynamical calculations 
of the collapse of molecular cloud cores with different initial rotation rates. From left to right, the different panels are for cloud cores with ^ = 0, 5 x 
10""^, 0.001, 0.005, 0.01, 0.04. The dashed lines give the evolution using the barotropic equation of state, while the solid lines give the evolution using 
radiation hydrodynamics. Time is set to zero at the end of the second dynamic collapse phase when the density reaches 10~^ g cm~^ which allows the length 
of the first hydrostatic core phases to be compared. It can be seen that the length of the first core phase increases with increasing rotation rate, but also that the 
use of radiation hydrodynamics and a realistic equation of state lengthens the first core phase relative to calculations that use the barotropic equation of state. 
Each calculation was performed with 10^ SPH particles, except the barotropic calculation with (3 = 0.04 which used 3 x 10^ particles. 



equation of state in all cases, by factors of 1.5 — 3 except for the 
most rapidly-rotating case. This is because the barotropic equation 
of state consistently underestimates the temperature of the gas, pro- 
viding less pressure support to the gas and, thus, allowing it to en- 
ter the second collapse phase earlier. The consistent temperature 
underestimate can be seen clearly in Fig. [6] which plots the maxi- 
mum temperature versus maximum density for each of the radiation 
hydrodynamical calculations and compares this to the barotropic 
equation of state (dashed black solid line). The lines from the radi- 
ation hydrodynamical calculations almost lie on top of one another, 
but from densities oflO~^^ to 10~^ g cm~^ the barotropic equa- 
tion of state underestimates the maximum temperature by about a 
factor of two. The only exception to this is the most rapidly-rotating 
calculation with /3 = 0.04. Here the maximum temperature at a 
given density is lower than in the other radiation hydrodynamical 
calculations and more similar to the barotropic equation of state. 
This is because in this case the first core is actually a torus and, 
therefore, the gas can cool more effectively. Returning to Fig.|5] we 
also see that it is the /3 = 0.04 case (rightmost panels) where the 
timescales between first and stellar core formation are most sim- 
ilar for the barotropic and radiation hydrodynamical calculations 
(differing only by about 15% rather than a factor of 1.5 — 3). 



One possible reason that the radiation hydrodynamical calcu- 
lations give greater temperatures during the first core phase than 
the barotropic calculations is that there may be considerable shock 
heating which is not taken into account with a barotropic equation 
of state. However, computing the entropy of the gas before the first 
core forms and in the bulk of the first core as it forms and evolves, 
we find it monotonically decreases in the radiation hydrodynami- 
cal calculations. The rate of decrease is rapid before the first core 
forms (when it is cooling rapidly and almost isothermal), but even 
after the first core forms the gas is loosing energy due to radiation. 
Only at the surface of the first core (around the accretion shock) 
does the entropy briefly increase before it radiatively cools. This is 
consistent with the recent calculations of Com mergon et aL] ( |2011| l 
who find that essentially all of the energy liberated in the accretion 
shock is radiated away. 

Instead, the reason the barotropic equation of state consis- 
tently underestimates the temperature in this density range is due to 
the approximation that 77 = 7/5 in this part of the evolution. In fact, 
molecular hydrogen (the dominant constituent) has a ratio of spe- 
cific heat capacities of 7 = 5/3 until it reaches ^ 100 K. Only then 
are the rotational degrees of freedom, which lower 7 to 7/5, excited. 
Again, this is apparent in Fig. [6] where it can be seen that the lines 
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Figure 6. The evolution of the maximum gas temperature versus maximum 
density for the radiation hydrodynamical collapse of molecular cloud cores 
with /3 = (soHd black line), 5 x 10~^ (dotted red curve), 0.001 (short- 
dashed green curve), 0.005 (long-dashed blue curve), 0.01 (dot-short dashed 
cyan curve), 0.04 (dot-long-dashed magenta curve). The barotropic equa- 
tion of state is given by the black short-dashed lines. The axisymmetric 
= 0) case is always the hottest at a given maximum density because 
the radiation is most effectively trapped. Conversely, the calculation that is 
the coolest has the highest rotation rate (/3 = 0.04). Each calculation was 
performed with 10^ SPH particles. 



from the radiation hydrodynamical calculations are steeper than the 
barotropic line in the temperature range ^ 20 — 150. Above this 
temperature, the lines are almost parallel (i.e. 7 ^ 7/5 for both), 
but the temperature offset (that originates in the 20 — 150 K range) 
persists. The barotropic equation of state could be improved for this 
part of the evolution by including a smooth transition from 7=1 
to 5/3 and then a transition from 5/3 to 7/5. However, this would 
still not capture all of the detail present in the radiation hydrody- 
namical equations (see below). 



3.2 Rotationally-unstable first hydrostatic cores 

If the first core is rotating rapidly enough that its own value of 
/3 > 0.274, the core is dynamically unstable to the growth of non- 
axisymmetric structure (|Bate|1998l|Machida et aL|2005l|Saigo &| 



Tomisaka'2006 , Saigo et al.'2008''Bate'2010';Machida et al.'2010! 
Tomida et al. 2010, Saigo & Tomisaka 2011, Machida & Mat- 
sumoto 2011). For the particular initial conditions used here, this 
occurs for the /3 = 0.001 — 0.01 cases (see Figs. [2] and [3]). The 
torus that forms in the /3 = 0.04 calculations is also dynamically 
unstable, but this case is somewhat different and will be discussed 
further below. 

In each of the /3 = 0.001 — 0.005 cases, the first core begins 
as an axisymmetric flattened pre- stellar disc, but after several ro- 
tations it develops a bar-mode (e.g. the second panel of the third 
row in each of Figures |2] and |3]). The ends of the bar subsequently 
lag behind and the bar winds up to produce a spiral structure. Spi- 
ral structure removes angular momentum from the inner parts of 
the first core via gravitational torques ( )Bate||1998| ), the effect of 
which can be seen in the evolution of density and temperature in 
Fig. [5] For example, in the /3 = 0.005 case, the slow increases 
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Figure 7. For the radiation hydrodynamical collapse of the molecular cloud 
core with /3 = 0.005, we give the evolution of the maximum gas tem- 
perature versus maximum density during the collapse (solid black curve) 
and snapshots of gas temperature versus density along the rotation (z) axis 
(long-dashed blue curves) and perpendicular to the rotation axis (dot-dashed 
red curves) at three different times. Note that the long-dashed blue lines at 
higher temperatures (later times) have a break in them. This is because there 
is not enough gas along the rotation axis at these densities to determine an 
accurate temperature. These densities occur within the accretion shock onto 
the first core. Although the evolution of the maximum temperature versus 
maximum density is similar to that given by the barotropic equation of state 
(black short-dashed lines), the gas temperature surrounding the protostar 
tends to be much hotter at a given density than predicted by the barotropic 
equation of state, particularly along the rotation axis where the tempera- 
tures can be more than an order of magnitude hotter. This calculation was 
performed with 3x10^ SPH particles. 



in central density and temperature after first core formation sud- 
denly accelerate with the onset of the spiral structure. This occurs 
at about 1300 yrs before stellar core formation for the radiation 
hydrodynamical calculation with P = 0.005 and about 600 yrs be- 
fore stellar core formation for the barotropic calculation. Similar 
accelerations are seen for both the /3 = 0.01 cases and also for the 
barotropic /3 = 0.001 case. Thus, the removal of angular momen- 
tum from the central regions of the core substantially accelerates 
the evolution of the first core towards the second collapse, which 
would take much longer to reach without the angular momentum 
redistribution (i.e. due to accretion and radiative cooling alone). 

In Figs.[2]and[3] the development of the spiral structure and the 
subsequent concentration of material towards the centre of the core 
due to the redistribution of angular momentum is clearly visible 
for the barotropic cases with /3 = 0.001 — 0.01 and the radiation 
hydrodynamical cases with /3 = 0.005 - 0.01. In Fig.|8] we also 
provide the density- weighted temperature, for comparison with the 
column density in Fig. [3] As gas is concentrated to the centre of the 
cores its temperature greatly increases. It is also apparent that the 
spiral shocks in the discs are hotter than the rest of the discs. 

The barotropic and radiation hydrodynamical evolutions are 
qualitatively similar in that both display the progression from an ax- 
isymmetric core, to bar instability, to a torus as the initial rotation 
rate of the molecular cloud core is increased. However, radiation 
hydrodynamical calculations are more resistant to the bar instabil- 
ity than the barotropic calculations. This is directly attributable to 
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Figure 8. Snapshots of the density- weighted temperature viewed parallel to the rotation axis during the evolution of the radiation hydrodynamical calculations 
of the collapse of molecular cloud cores with different initial rotation rates. From top to bottom, the different rows are for cloud cores with /3 = 5 x 
IQ-^, 0.001, 0.005, 0.01, 0.04. Note that the spatial scale is different for each row, with each panel measuring 1600^^^^ AU across (i.e. from 36 to 320 AU). 
The temperature sale differs between the top three rows and the lower two rows in order to enhance the low-temperature structure in the cases with greater 
rotation. The free-fall time of the initial cloud core, tff = 1.8 x 10-*^^ s (56,500 yrs). Each calculation was performed with 10^ SPH particles, except for the 
/3 = 0.001 and (3 = 0.005 cases which used 3 x 10^ SPH particles. The evolution is almost identical when using 10^ or 3 x 10^ particles, but the latter are 
slightly more detailed (see Appendix A). 
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the difference in the thermal evolution of the gas that was discussed 
above. Since the gas is somewhat hotter in the radiation hydrody- 
namical calculations, the first core is larger (more 'puffy') and has 
a lower value of /3 for a given /3 value of the initial molecular cloud 
core. |Tomida et al.| { |2010| ) also noticed that the first cores in their ra- 
diation magnetohydrodynamical calculations had higher entropies 
and larger sizes than when they used a barotropic equation of state. 
It should be noted, however, that with radiation hydrodynamics 
there is no one-to-one relation between temperature and density 
(see also |Boss et al.|2000[|Whitehouse & Bate|2006[rTomida et al.| 
|2010| l. This is illustrated in Fig.|7] where we plot temperature ver- 
sus density from various snapshots from the /3 = 0.005 calculation 
and compare this to both the run of maximum temperature versus 
maximum density during the calculation and the barotropic equa- 
tion of state. The red dot-dashed lines are for cuts perpendicular 
to the rotation axis along the x-axis (i.e. in the plane of the disc), 
while the blue long-dashed lines display the temperature along the 
rotation axis. Not only is the run of maximum temperature versus 
maximum density always greater than or equal to the temperature 
from the barotropic equation of state, but the gas both in the mid- 
plane of the disc and along the rotation axis at the same density is 
even hotter. This is particularly true of the gas along the rotation 
axis, which is heated by the accretion shock at the surface of the 
first core. |Tomida et al.| ( |201 0) give a very similar plot for a snap- 
shot from one of their simulations and find very similar behaviour. 

These hotter temperatures in the radiation hydrodynamical 
calculations mean that while the transition from rotationally sta- 
ble to dynamically unstable first core occurs in the range /3 = 
using the barotropic equation of state, the transi- 
tion occurs at /3 ^ 0.001 with radiation hydrodynamics (this case 
undergoes an extremely weak instability, just visible in the fourth 
panel of the second row of Fig.|3]). Similarly, while the barotropic 
calculation with /3 = 0.01 manages to produce a torus-shaped first 
core, the radiation hydrodynamical calculation with (3 — 0.01 is 
still definitely disc shaped rather than torus shaped. 

As mentioned above, for the highest initial rotation rates (/3 = 
0.04) the first core is actually a torus or ring-like structure (e.g. 
INorman & Wilson|1978||Cha & Whitworth|20031 [Machida et al.| 
'2005). Such a configuration is highly unstable to non-axisymmetric 
perturbations and, indeed, as is clearly visible in Figs. [2] and [3] the 
rings rapidly fragment into four objects. Such a configuration is 
highly chaotic, and symmetry is broken quickly (due to truncation 
error and the use of a tree-structure to calculate gravity in the cal- 
culations). In the radiation hydrodynamical calculation, two of the 
fragments merge to produce a triple system, while in the barotropic 
calculation all four fragments survive (at least until the calculation 
was stopped). Each of the fragments follows its own evolution to- 
ward the second collapse phase and stellar core formation. The cal- 
culations were stopped soon after the first fragment in each calcu- 
lation produced a stellar core. 

Finally, as mentioned above, first cores may evolve into pre- 
stellar discs with radii ranging from ^ 5 to ^ 100 AU before a stel- 
lar core forms (Fig.jsj). Those with radii greater than ^ 10 AU are 
produced due to the angular momentum transport that occurs dur- 
ing the dynamical rotational instability. A further consequence is 
that all star formation should go through at least a brief phase when 
the disc mass is greater than the stellar mass. Such discs are grav- 
itationally unstable and may evolve through spiral density waves 
(e.g. the third row of Figsjsjand [8]) and/or fragmentation (e.g. the 
fourth rows of Fig.|3]and [s])? This is discussed further in Section|4] 

In conclusion, in switching from the simple barotropic equa- 
tion of state to a realistic equation of state with radiation hydrody- 



namics, the qualitative evolution of the first core and its dependence 
on the initial rotation rate of the molecular cloud core is identi- 
cal. However, the gas temperature in the more realistic calculations 
tends to be hotter and is not only a function of density: for example, 
the gas is significantly hotter in the accretion shock at the surface 
of the first cores. Quantitatively, the higher temperatures mean that 
the critical values of the initial molecular cloud core rotation rates 
required for bar instability of the first core, or the transition to a 
torus geometry is somewhat higher. The value of /3 must be 50% 
greater with radiation hydrodynamics, which translates into a rota- 
tion rate which is ~ 25% greater (since ^ (x , where is the 
angular frequency of the molecular cloud core). 

3.3 The effect of stellar core formation 

Although the evolution of the first core or pre- stellar disc is qualita- 
tively the same when computed with a barotropic equation of state 
or radiation hydrodynamics, as |Bate| ( [2010| ) showed the evolution 
subsequent to the formation of the stellar core is qualitatively dif- 
ferent. When using a barotropic equation of state, the formation of 
the stellar core deep within the optically-thick disc has no effect on 
the temperature of the gas further out in the disc because its temper- 
ature is set purely according to its density. However, with radiation 
hydrodynamics, the situation is completely different. 

When the second collapse occurs and produces the stel- 
lar core, the gravitational potential energy that is released is ^ 
GMsc/Rsc = 4 X 10^^ erg. Since the stellar core is in virial equi- 
librium, approximately half of this energy is radiated away. More- 
over, the stellar core rapidly begins to accrete from the first core, 
reaching a mass of 6 Mj in only a few years (i.e. more quickly 
than the dynamical timescale of the large-scale disc; see bel ow). 
This increases the total energy released to 3 x 10^^ erg. Bate] 
( |2010| ) compared this energy with the binding energy of the disc in 
the /3 = 0.005 calculation. Just before the onset of the second col- 
lapse, the binding energy of the pre-stellar disc is only 4 x 10^^ erg 
(estimated as GM^/Rd with Md ^ 0.18 M© and taking a 
mean 'spherical' radius of Rd ~ 15 AU). Thus, when the second 
collapse occurs, the disc suddenly finds itself irradiated from the 
inside by an energy source emitting a substantial fraction of the 
binding energy of the disc itself. Because the disc is extremely op- 
tically thick, this energy is temporarily trapped in the centre of the 
disc and heats the gas dramatically, sending a weak shock wave out 
along the midplane of the disc at a speed of a few km s~^. How- 
ever, perpendicular to the disc, the effect is even more dramatic. 
Because there is less material along the rotation axis, the hot gas 
finds it easiest to break out in this direction and a bipolar outflow 
is launched. Whereas the wave within the disc decays as it travels 
leaving the bulk of the disc gravitationally bound, the gas form- 
ing the bipolar outflow has velocities up to 10 km s~^ and travels 
out into the infalling envelope to distances in excess of 50 AU in 
less than 50 years. Using two-dimensional radiation hydrodynam- 
ical calculations, ^Schonke & Tscharnuter (2011 ) have recently re- 
ported similar behaviour, although the outflow in their simulations 
is less bipolar. By disguarding the central regions of one of their 
calculations, they were able to follow the outflow for a few tens of 
thousands of years and found that after reaching ^ 500 AU, the 
material fell back in to reform a disc. 

| Bate| ( [207^ mainly discussed the /3 — 0.005 calculation as 
typical example. Here we examine how the effect of stellar core 
formation on the surrounding disc depends on the initial rotation 
rate of the molecular cloud core and, thus, on the degree of rotation 
that the first core/pre-stellar disc has. Figs. |9] to [14] illustrate the 
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Figure 9. Snapshots of the column density viewed parallel to the rotation axis during the evolution of the radiation hydrodynamical calculations of the collapse 
of molecular cloud cores with different initial rotation rates. The Shockwaves propagating outwards through the discs following stellar core formation are 
clearly visible. From top to bottom, the different rows are for cloud cores with /5 = 5 X 10-4, 0.001, 0.005, 0.01, 0.04. Note that the spatial scale is different 
for each row, with each panel measuring IGOOy^ AU across (i.e. from 36 to 320 AU). The free-fall time of the initial cloud core, tff = 1.8 x 10-*^^ s 
(56,500 yrs). Each calculation was performed with 10^ SPH particles, except for the (3 = 0.001 and (3 = 0.005 cases which used 3 x 10^ SPH particles. 
The evolution is almost identical when using 10^ or 3 x 10^ particles, but the latter are slightly more detailed (see Appendix A). 



evolution of calculations with /3 = 5 x 10~^ — 0.01 following 
stellar core formation. The /3 = case is not discussed further since 
the situation is as |Larson| ( |1969| ) described it. This case remains 
spherically- symmetric, and although the stellar core irradiates the 
first core from within, the radiation liberated is not sufficient to stop 
the spherically- symmetric accretion flow onto the stellar core. We 
do not discuss the /3 = 0.04 case following stellar core formation 
since each of the three fragments is qualitatively similar to the first 
core obtained in the (3 — 0.001 case and, therefore, we assume that 
each of these cores would evolve in a similar manner. 



Figs. [9] and [T0| provide snapshots of the column density in 
the directions parallel to the rotation axis and perpendicular to the 
rotation axis, respectively. The Shockwave propagating outwards 
through each of the first cores/discs is clearly visible in Fig.|9] For 
the /3 = 5 X 10"^ and 0.001 cases, the Shockwave reaches the 
outer edge of the disc before the calculations are stopped. For the 
P = 0.005 and 0.01 cases, we were not able to follow the calcula- 
tions this long and the Shockwave was only followed 10 — 15 AU. In 
Fig. [To| both the propagation of the Shockwaves through the discs 
and the biploar outflows launched perpendicular to the discs are 
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Figure 10. Snapshots of the column density viewed perpendicular to the rotation axis during the evolution of the radiation hydrodynamical calculations of the 
collapse of molecular cloud cores with different initial rotation rates. The bipolar outflows, launched after stellar core formation, are clearly visible. From top 
to bottom, the different rows are for cloud cores with /5 = 5 x 10~^, 0.001, 0.005, 0.01, 0.04. Note that the spatial scale is different for each row, with each 
panel measuring 1600^^ AU across (i.e. 36, 50, 114, 160, or 320 AU). The free-fall time of the initial cloud core, tff = 1.8 x 10^^ ^ (56,500 yrs). Each 
calculation was performed with 10^ SPH particles, except for the /3 = 0.001 and ^ = 0.005 cases which used 3x10^ SPH particles. The evolution is almost 
identical when using 10^ or 3 x 10^ particles, but the latter are slightly more detailed (see Appendix A). 



clearly visible, the latter reaching to distances of 15-35 AU. The 
furtherest an outflow was followed was to a distance of 60 AU in 
the P — 0.005 case, approximately 50 years after the outflow began 
( |Bate|2010| ). 

Each of Figs. [TT] to [13] give density, velocity, and temperature 
profiles, both in the disc plane and perpendicular to the disc plane 
(i.e. along the rotation axis), at four characteristic times during the 
evolution following stellar core formation. They also provide the 
radial mass profiles. Fig. [TT] shows the evolution of the barotropic 
P = 0.005 calculation which can be compared with the results 



from the radiation hydrodynamical calculations in Figs. [T2|to[T3| 
with /3 = 5xl0~^ and 0.005. We do not provide figures for /3 = 
0.001 and 0.01 since they are qualitatively similar to the cases with 
/3 = 5 x 10"^ and 0.005, respectively. 

In all calculations, the stellar core forms with a radius of 
Rsc ^ 0.01 AU (i.e. 2i?0) and a mass of a few Mj (e.g. Fig. [TT] 
top-left and top-centre panels). Using the barotropic equation of 
state (Fig.[TT]), the radii of the first core (actually a disc) and stellar 
core are clearly visible in the top-centre panel which gives the ra- 
dial velocity profiles just after stellar core formation. Gas falling on 



14 M.R. Bate 




-2 2 

Radius, logjo r [AU] 




-2 2 

Radius, logj^ r [AU] 



-2 2 

Radius, logj^ r [AU] 



-2 2 

Radius, log,(j r [AU] 





_ , 1 , , , 1 , , 


1 ' ' ' - 










- , 1 ,v , , 1 , , 


1 , , , - 



Radius, log,Q r [AU] 



2 
Radius, logj^ r [AU] 



-2 2 

Radius, logj^ r [AU] 



-2 2 

Radius, log,o r [AU] 







-2 2 

Radius, log,^ r [AU] 




-2 2 

Radius, log,(j r [AU] 




-2 2 

Radius, log,(j r [AU] 



Figure 11. The evolution of the barotropic /5 = 0.005 calculation of the collapse of a molecular cloud core following stellar core formation. Each of the four 
rows shows the state of the protostar at different time. The left panels provide the radial density profile perpendicular to the rotation axis (solid line, averaged 
in azimuth), the density profile along the rotation axis (dashed line), and the cumulative mass profile (dotted line). The centre panels give the radial velocity 
profiles perpendicular to (solid line) or along (dashed line) the rotation axis. The right panels show the radial temperature profiles perpendicular to (solid line) 
or along (dashed line) the rotation axis. Breaks in the lines occur when there is not sufficient gas resolution to define an accurate measurement. 



to the first core from the envelope at speeds of 2 — 4 km s~^ is 
visible at radii of ^ 10 — 30 AU, while gas falling onto the stellar 
core much more rapidly (up to 7 — 10 km s~^) is apparent at radii 
of ^ 0.01 — 0.3 AU. In less than a year (second and third rows 
of Fig. [TT]), an inner disc builds up around the stellar core (radius 
^0.07 AU in the centre panel of the second row, and 0.2 AU in 
the centre panel of the third row) and the radial extent of the region 
undergoing dissociation and collapse on to the stellar core and inner 
disc is reduced. Several years after the stellar core has formed (bot- 
tom row of Fig.[TTJ, the outer edge of the inner disc has completely 
merged into the outer disc (the remnant of the first core) and the 
azimuthally- averaged radial velocity in the disc plane is approxi- 
mately zero out to ~ 30 AU. The formation of an inner disc around 
the stellar core, which grows in radius as gas with higher angular 
momentum falls in through the region of molecular hydrogen dis- 
sociation, and the eventual merger of the inner and outer discs is 
just as was first discussed by |Bate|fl998| ), and more recently stud- 



ied by |]V[achida et aLlpOlOT ) and [IMachida & ]V[atsumotol { |2011| ) 
who also included magnetic fields. 

Using radiation hydrodynamics, the evolution is initially sim- 
ilar to the barotropic cases. For example, take the /3 = 5 x 10~^ 
case (Fig. [12]). In the top-centre panel, the first core and stellar cores 
are clearly visible in the radial velocity profiles due to the accretion 
shocks at their surfaces. An inner disc builds up around the stellar 
core as material with higher angular momentum falls inwards and 
the radial extent of the dissociating, collapsing region decreases 
(centre panel of the second-row). 

However, as the inner disc grows in radius and begins to merge 
with remnant of the first core, a bipolar outflow is launched verti- 
cally at about 5 km s~^ (dashed line in the centre panel of the third- 
row of Fig. [12]), while a lower- velocity shock propagates outward 
along the disc mid-plane at 1 km s~^ . All of this evolution occurs 
within just 4 — 5 years of the stellar core forming. At 20 years after 
stellar core formation, the outflow has travelled more than 20 AU 
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Figure 12. The evolution of the radiation hydrodynamical (3 = 5xl0~^ calculation of the collapse of a molecular cloud core following stellar core formation. 
Each of the four rows shows the state of the protostar at different time. The left panels provide the radial density profile perpendicular to the rotation axis (solid 
line, averaged in azimuth), the density profile along the rotation axis (dashed Hne), and the cumulative mass profile (dotted Hne). The centre panels give the 
radial velocity profiles perpendicular to (solid line) or along (dashed line) the rotation axis. The right panels show the radial temperature profiles perpendicular 
to (solid line) or along (dashed line) the rotation axis. 



and the Shockwave in the disc has burst out of the edge of the disc 
(centre panel of the bottom-row of Fig. [12} see also Fig.[TO|. 

The outflow and the Shockwave do not vary greatly with the 
degree of rotation of the first core/pre- stellar disc (c.f. Figs, and 
[13]). In all cases, the outflow is launched (vertically) at a radius of 
^ 0.5 AU from the stellar core, just before the outer edge of the 
inner disc merges (in the disc plane) with the remnant of the first 
core/pre- stellar disc at ^ 1 AU and the infalling region of disso- 
ciating molecular hydrogen disappears (e.g. the centre panel of the 
second row in Fig. [TsJ. The peak speed of the outflow is slightly 
larger when the first core is rotating more rapidly, ranging from 
^ 5 km s~^ for the /3 = 5 x 10~^ case to ^ 10 km s~^ for the 
/3 = 0.01 case. The outflows have been followed for long enough 
in the /3 = 0.001 and /3 = 0.005 calculations (to distances of 30 
and 60 AU, respectively) that speeds of the outflows were clearly 
seen to decay (c.f. the dashed lines in the centre panels of the third 
and fourth rows of Fig. [13]). This outflowing material will not leave 
the system; it will eventually stall and fall back onto the proto- 



star and disc { [Schonke & Tscharnuter|20lT] ). It is worth noting that 
because the flux-limited diffusion method used to calculate the ra- 
diative transfer in the calculations presented here assumes that the 
material radiates with maximal efficiency, it may overestimate the 
cooling and hence underestimate the temperatures in the optically- 
thin outflowing material. Thus, the outflow may be even stronger 
with a more accurate method of radiative transfer. 

The Shockwave propagating outwards through the disc travels 
with a speed of 1 — 2 km s~^ in all cases. In the /3 = 5 x 10~^ and 
/3 = 0.001 calculations, the Shockwave is followed until after it has 
reached the edge of the first core/pre- stellar disc. The combination 
of the outflow and the outward propagating wave in the disc plane 
evacuate most of the gas from within a few AU of the protostar, 
leaving the stellar core surrounded by a broad torus, rather than a 
disc. This is apparent in the column density plots in Fig. [9] as well 
as the mass and density profiles in the bottom rows of Figs. [T2] to 
[T3] These holes would eventually fill in again as the disc evolves. 
In the low-resolution calculations, the mass resolution is poor and 
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Figure 13. The evolution of the radiation hydrodynamical (3 = 0.005 calculation of the collapse of a molecular cloud core following stellar core formation. 
Each of the four rows shows the state of the protostar at different time. The left panels provide the radial density profile perpendicular to the rotation axis (solid 
line, averaged in azimuth), the density profile along the rotation axis (dashed fine), and the cumulative mass profile (dotted fine). The centre panels give the 
radial velocity profiles perpendicular to (solid line) or along (dashed line) the rotation axis. The right panels show the radial temperature profiles perpendicular 
to (solid line) or along (dashed line) the rotation axis. 



the 'holes' around the stellar cores are in fact holes (devoid of SPH 
particles). However, in the high resolution /3 = 0.001 and (3 — 
0.005 cases (which each use 3 x 10^ particles), it is found that 
a small amount of gas does remain within the holes, and this gas 
continues to slowly accrete onto the stellar core. 



In Fig. 



14 



we plot the mass with density > 10~ g cm~ 
versus time for the /3 = 0.001 and /3 = 0.005 cases with dif- 
ferent resolutions. It can be seen that the mass of the stellar core 
at which the feedback dramatically curtails the accretion decreases 
from 10 — 11 down to 6 Mj with increased resolution. For a few 
years, the stellar cores grow at a rate of ^ 10~^ M© yr~^ (even 
with the highest resolutions). With low-resolution (^ 1 x 10^ SPH 
particles), the accretion onto the stellar cores ceases entirely after 
the launching of the outflows due to the formation of the 'holes' in 
the inner few AU of the discs. However, using 3 x 10^ SPH parti- 
cles, it is found that the stellar cores do continue to accrete during 
this time, but at much reduced accretion rates of 2 x 10~^ M© yr~^ 
for the /3 = 0.001 calculation and 1 x 10~^ M© yr"^ for the 



(3 — 0.005 calculation (measured between 20 and 50 years after 
stellar core formation; Fig. [14]). 



Although the rate of convergence with increasing resolution is 
relatively slow, strong outflows are launched in all cases (see Ap- 
pendix A). The slow convergence is due to the strong interplay be- 
tween the energy released by the formation of the stellar core, and 
the launching of the outflow which decreases accretion and, thus, 
reduces the energy that is released. If the resolution is poor, the 
stellar core accretes more material before the energy released feeds 
back and manages to stop the accretion, whereas with higher res- 
olution, the accretion of a small amount of mass can feed thermal 
energy into the infalling material more quickly and, thus, inhibit 
further accretion. 
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Figure 14. The mass of the stellar core (gas with density > 10 ~ g cm~^) measured from the time of its formation. The two panels shows the results for 
the radiation hydrodynamical (3 = 0.001 (left) and (3 = 0.005 (right) molecular cloud cores performed using 10^ (long-dashed line), 3 x 10^ (short-dashed 
line), 10^ (dotted line), and 3 x 10^ (solid line) particles. With low resolutions the radiative feedback stops accretion onto the stellar cores, but for the highest 
resolutions the accretion continues at a low level (1 — 2 x 10~^ Mq yr~^). 



4 DISCUSSION 

4.1 The lifetime of the first core/pre-stellar disc phase 

As mentioned in Section [3l2l the introduction of radiative transfer 
and a realistic equation of state increases the lifetimes of the first 
core phase compared to that obtained using our barotropic equa- 
tion of state by factors of 1.5-3 (Fig.|5]). This is important, because 
the longer the first core phase lasts, the more easy it will be to ob- 
serve. The lifetimes obtained using radiative transfer range from 
^ 400 years (with no rotation) to ^ 3000 years (for /3 = 0.01), 
whereas using the barotropic equation of state they ranged from 
^ 100 to ^ 1500 years. As discussed above, the difference is due 
to the higher temperatures (and thus higher pressures) that are ob- 
tained using the realistic physics rather than the barotropic equa- 
tion of state, which slows the evolution towards the second col- 
lapse phase. This lengthening of the lifetimes with radiation hydro- 
dynamics compared to barotropic calculations was also seen by To- 
|mida et aL|p010| ). The lengthening occurs for both the non-rotating 
first cores, and in the cores that undergo non-axisymmetric instabil- 
ities. In the latter, once the first cores or pre- stellar discs develop the 
bar instability they very quickly evolve to high central densities and 
temperatures and undergo collapse to stellar densities when com- 
puted using the barotropic equation of state (Fig. [5] 4th and 5th 
panels across where the maximum density increases rapidly from 
10~^^ to 10~^ g cm~^, and the maximum temperature from 100 
to 2000 K). Using the more realistic physics, this evolution takes 
3-8 times longer. This is because the hotter temperatures weaken 
the strength of the spiral arms generated by the instability (c.f. the 
third rows of Figs. [2] and |3]) and, thus, decrease the rate of angular 
momentum transport within the pre- stellar disc. The removal of the 
angular momentum from the central regions of the disc removes 
rotational support from the gas (B ate,1998j ), allowing it to contract 
and heat up much more rapidly than it would due to accretion alone, 
eventually triggering second collapse. 

The pre- stellar disc phase will be more easily observed in 
more rapidly-rotating molecular cloud cores because the lifetime 
of the phase is longer. Resolved observations of the pre- stellar disc 
phase (e.g. using ALMA) would also be easier and more interesting 
for the more rapidly-rotating cases because the disc will be much 



larger (up to 50 — 100 AU in radius rather than 5 — 10 AU 
for rotationally- stable first cores) and substructure may be visi- 
ble (i.e. the presence of spiral density waves). Recently, Sai go &| 
|Tomisaka|(2011 ) investigated the observability of pre-stellar discs 
using ALMA and concluded that ALMA should be able to image 
objects in the nearest star-forming regions. For the particular initial 
conditions used in this paper, the lifetimes of as a function of the 
initial cloud rotation rate are plotted in Fig. [15] As pointed out re- 
cently by |Tomida et al.|p010| ), the lifetime of the first core phase 
also depends on the mass of the molecular cloud core such that the 
lifetimes are longer in lower-mass molecular cloud cores. For ini- 
tial core masses of^ 0.1 M©, the accretion onto the first core/pre- 
stellar disc is not sufficient to drive the object into the second col- 
lapse phase. Instead, the central regions of the first core/disc evolve 
to higher temperatures primarily through radiative cooling. As the 
first core radiates energy, it contracts and heats up, eventually ex- 
ceeding 2000 K when the dissociation of molecular hydrogen be- 
gins, but this process can take in excess of 10^ years. Thus, the first 
core/pre-stellar disc phase is more likely to be observed for molec- 
ular cloud cores that have higher rotation rates and lower masses. 

Despite this, observing the pre-stellar disc phase with ALMA 
will be challenging. Lifetimes of even 4 x 10^ — 10^ years are 
still only 2-5% of the estimates of Class lifetimes of 2 x 10^ yrs 
( |Hatchell et al.|2007| [Evans et aL|2009| l, so that for every 100 Class 
objects there should only be a few pre-stellar discs. The first place 
to look is likely to be the so-called very low luminosity objects 
(VeLLOs) with luminosities of less than 0.1 Lq (e.g. | Young et aT] 
[20041 ICrapsi et al.||2005[ [Huard et al.||2006l |Bourke et al.||2006| ), 
some of which have recently been proposed as candidate first cores 
( |Chen et al., 20 10, Enoch et al.|2010j . 



4.2 Disc to stellar mass ratios 

As first discovered by |Bate| ( |1998| ), a rotating first core actually 
leads to a disc that forms before the star. These discs can range from 
^ 5 AU in radius for very slowly rotating first cores (c.f. |Larson| 
7969t to ^ 100 AU for cores that undergo rotational instabilities 
UMachida et aL|2010|[Batei20T0|[Machida & Matsumoto|2011| ). 
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Because the disc forms first, there is always a phase of star for- 
mation when the disc to stellar mass ratio is greater than unity. For 
example, in the original calculation of Bate ( 1998 ), the disc had a 
mass of 0.08 Mq when the stellar core formed with an initial mass 
of only 1.5 Mj, giving a disc to stellar mass ratio ^ 50. Although 
some quite massive discs have been discovered around some young 
stars (e.g. the 0.8 M© disc around a 3.5 Mq star; |Hamidouche| 
|20T ) such cases are rare (e.g. Boissier et al. 2011) and the usual 
assumption is that the circumstellar disc mass is always likely to be 
less than the stellar mass. While this is almost certainly true later 
in the evolution of a protostar (e.g. Class I or Class II low-mass ob- 
jects and Herbig Ae/Be stars), in the very earliest phases the disc to 
star mass ratio can greatly exceed unity. 

With the initial conditions used in this paper, the disc mass 
when the stellar core forms increases with the rotation rate of the 
progenitor molecular cloud core. This is because a greater fraction 
of the 1 M0 molecular cloud core has sufficient angular momentum 
to settle outside the typical 5 AU radius of a non-rotating first core. 
Because more mass settles away from the centre of the object, the 
accretion from the infalling envelope is less effective at increasing 
the central temperature (i.e. mass receives support from rotation as 
well as thermal pressure). Thus, the first core phase lasts longer and 
the masses of the first cores/pre- stellar discs at the time of stellar 
core formation range from ^ 0.03 Mq for /3 = to 0.22 Mq 
for /3 = 0.01 (see Fig. [15] and the dotted lines in the left panels of 
Figs. [12] to [13]). In the most extreme case we have modelled here, 
when the stellar core forms with an initial mass of ^ 2 Mj, the 
disc is more than 100 times more massive than the star! Note that 
the initial rotation rate of the molecular cloud core is not extreme 
for this case. [Goodman et al.| ( |1993| ) find that the values of /3 for 
observed molecular cloud cores range from 0.002 to > 1, with typ- 
ical values of /3 ^ 0.02. Therefore, it is quite possible that when 
pre-stellar discs are observed (e.g. with ALMA), they could have 
substantial masses (e.g. up to 1/4 of the mass of an entire l-M© 
molecular cloud core) without a stellar core having formed. 



4.3 Disc fragmentation 

The potential for high disc to stellar mass ratios also implies that 
some of these discs may fragment very early in the star forma- 
tion process. Indeed, for our /3 = 0.01 case (and, of course, the 
/3 = 0.04 case), this is exactly what happens, with the disc frag- 
menting to produce two additional 'first cores' even before the orig- 
inal first core/disc has undergone second collapse to produce a stel- 
lar core. The fate of such fragments is of course not clear. Multiple 
fragments may merge with each other before they each undergo 
second collapse. They may also spiral in to the central object and 
be disrupted before they undergo second collapse. However, those 
that survive will produce the seeds of binary or multiple stellar sys- 
tems. 

Whitehouse & Bate I pOPg]), [Krumholz 



As pointed out by 



( [2006] ), [Krumholz et al.| p007| ), |Bate| (|2009) and "Offne r eTal 
P009|), radiative feedback from newly formed proto stars heats their 



surrounding gas and inhibits fragmentation. In star cluster forma- 
tion simulations (e.g. [Bate|2009l [Offner et al.|2009j ) this radiative 
heating dramatically reduces the number of objects formed com- 
pared to using a barotropic equation of state. Some of this reduc- 
tion comes from a reduction in the frequency of disc fragmentation 
because the discs are heated and stabilised by the radiation from 
their central objects ( |Bate|2009| ). 

However, as shown here in the P = 0.01 case, fragmentation 
of the pre-stellar disc resulting from a rapidly-rotating first core 
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Figure 15. The dependence of the pre-stellar disc (first core) properties 
on the initial rotation rate of the progenitor molecular cloud core, /5. We 
plot the masses and radii of the pre-stellar cores, measured just before the 
second collapse to form the stellar core occurs, and the lifetime of the pre- 
stellar disc (i.e. the time between the formation of a first hydrostatic core 
and the formation of the stellar core). The mass, radius, and lifetime all 
increase with larger initial molecular cloud core rotation rates. In the more 
rapidly-rotating cases, the pre-stellar disc accumulates a significant fraction 
of the total cloud mass before the stellar core forms. Note that most rapidly- 
rotating /5 = 0.04 case has been excluded because it forms four fragments 
rather than a single pre-stellar disc. 



phase can occur before the central stellar core is formed. Thus, even 
if radiative feedback can completely stop disc fragmentation after 
the stellar core forms, disc fragmentation before stellar core for- 
mation may be a common route to forming binary and multiple 
systems. 



4.4 Implications of radiative feedback following stellar core 
formation 

The radiative impact of stellar core formation on the surrounding 
disc, driving a shock wave through the disc and launching an out- 
flow may have several wider implications. These were discussed 
by [Bate|pOTO| , so we only briefly list them here. First, in future 
studies it will be important to examine the relative roles of radia- 
tive heating due to stellar core formation and magnetic fields in the 
launching of a protostellar jet. Second, although the outburst and 
the dramatic decrease in the stellar core accretion rate discussed in 
this paper is a transient associated with stellar core formation, it 
might reoccur if the accretion rate onto the stellar core increases 
again and, thus, may potentially be a source of episodic accretion 
and outbursts. Third, such outbursts may contribute to the produc- 
tion of the chondrules found in meteorites ( [Grossman et al.|1988] ). 



5 CONCLUSIONS 

We have presented results from three-dimensional radiation hy- 
drodynamical calculations that follow the collapse of a molecular 
cloud core through the formation and evolution of the first hydro- 
static core and beyond the formation of the stellar core. We find 
the evolution before the formation of the stellar core is qualita- 
tively similar to that found in the past using barotropic equations of 
state. As found in past studies, the evolution of the first hydrostatic 
core depends on the initial rotation rate of the cloud. Non-rotating 
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or slowly rotating clouds produce small spherical or weakly flat- 
tened first cores with radii of ~ 5 — 10 AU. More rapidly ro- 
tating clouds (/3 ^ 0.001 — 0.01) produce highly-oblate rapidly- 
rotating first cores that undergo dynamic rotational instabilities to 
produce large discs (radii ^ 10 — 100 AU) before the formation 
of a stellar core. These objects are better described as pre- stellar 
discs than first cores. Yet more rotation typically results in frag- 
mentation, via a torus shaped first core if the initial cloud is ax- 
isymmetric. Quantitatively, using radiation hydrodynamics with a 
realistic equation of state produces first cores that are somewhat 
hotter, slightly more stable to rotational instabilities, and longer- 
lived (by factors of 1.5 — 3) than using our barotropic equation of 
state. 

The masses and radii of the pre- stellar discs produced before 
stellar core increase with the initial rotation rate of the molecular 
cloud core. Their lifetimes also increase with rotation rate, and for 
lower-mass molecular cloud cores (see also [Tomida et al.||2010] l. 
With high rotation rates, these pre-stellar discs can have radii up to 
^ 100 AU and contain in excess of 0.2 M© of gas (up to 1/4 of 
the mass of the entire molecular cloud core) and exist for several 
thousands of years before stellar core formation. Such objects may 
be resolvable by ALMA. Since initial mass of the stellar core which 
forms within these discs is just a few Jupiter-masses, the disc-to- 
star mass ratios can exceed 100 for a short period of time. 

Fragmentation of a pre-stellar disc prior to stellar core forma- 
tion may be an effective way to produce binary and multiple star 
systems because the radiative feedback from the stellar core which 
may inhibit later fragmentation of a massive disc is absence. 

After the formation of a stellar core within the disc, dissociat- 
ing gas from the inner regions of the disc falls towards the stellar 
core and builds a smaller disc around the stellar core. This inner 
disc increases with radius as material with higher angular momen- 
tum falls in, and eventually merges with the inner edge of the outer 
disc. 

Although this inner/outer disc structure and the earlier evo- 
lution of the first core phase is similar for barotropic and radiation 
hydrodynamical calculations, we find that the subsequent evolution 
is qualitatively different in radiation hydrodynamical calculations. 
In barotropic calculations, the formation of the stellar core deep 
inside the first core (or pre-stellar disc) has no effect on the sur- 
rounding disc because the temperature of the gas is simply set by 
the density of the gas. However, with radiation hydrodynamics, the 
energy released by the formation of the stellar core within the op- 
tically thick disc is similar to the binding energy of the disc. This 
heats region surrounding the stellar core and, at approximately the 
same time as the inner disc and outer discs merge, a shock wave 
forms in the merger region and propagates outwards through the 
disc, while a bipolar outflow is launched perpendicular to the rota- 
tion axis, dramatically decreasing the accretion rate on to the stellar 
core. The properties of the outflow do not seem to vary greatly with 
the initial rotation rate of the molecular cloud core or the pre-stellar 
disc, but the speed of the outflow is somewhat greater for higher ro- 
tation rates and the Shockwave which propagates outward through 
the disc reaches the outer edge of the disc more quickly for smaller 
discs. It remains to be seen how the inclusion of magnetic fields 
alters this thermally-driven outburst. 
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APPENDIX A: NUMERICAL CONVERGENCE OF THE 
RADIATION HYDRODYNAMICAL CALCULATIONS 



As mentioned in Section |2.2| for the different initial conditions 
and equations of state discussed in this paper, most of the calcula- 
tions were performed with several different numerical resolutions 
to check for numerical convergence of the results. The numbers of 
SRH particles used were 1 x 10^ 3 x 10^ 1 x 10^ and 3 x 10^ 
particles, although as summarised in Table[T]not all cases were per- 
formed with all resolutions. The only type of calculation to be per- 
formed with all four resolutions was the /3 = 0.005 case with ra- 
diation hydrodynamics, but the evolution of three other initial con- 
ditions performed with radiation hydrodynamics were studied with 
two (P = 0.01 and 0.04) or three (/3 = 0.001) different resolutions. 
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Figure 16. Examining the numerical convergence of the lifetime and evo- 
lution of the first core/disc. We plot the time evolution of the maximum 
density (upper panels) and gas temperature (lower panels) during the radi- 
ation hydrodynamical calculations of the collapse of molecular cloud cores 
with different initial rotation rates. From left to right, the different panels are 
for cloud cores with (3 = 0.001, 0.005, 0.01, 0.04. The different line types 
are for calculations performed using 1 x 10^ (long-dashed line), 3 x 10^ 
(dotted lines), 1 x 10^ (short-dashed lines), and 3 x 10^ (solid lines) SPH 
particles. Time is set to zero at the end of the second dynamic collapse phase 
when the density reaches 10~^ g cm~^ which allows the length of the first 
hydrostatic core phases to be compared. It can be seen that the lifetime of 
the first core phase is well converged for all except the lowest resolution cal- 
culation for the ^ = 0.005 case, and the /5 = 0.001 case. This is because 
the lifetime of the first core/disc is determined by the dynamic rotational in- 
stability for /5 ^ 0.005, whereas for /3 = 0.001 viscous evolution probably 
plays a significant role. 



In this appendix, we discuss how the results depend on numerical 
resolution. 

As numerical resolution was increased, all of the clouds took 
slightly longer to collapse. However, more important for the re- 
sults presented in the main sections of this paper is how the life- 
time and evolution of the first core or pre- stellar disc phase de- 
pends on numerical resolution. In Figure [16] we plot the evolution 
of maximum density and temperature versus the time before stel- 
lar core formation (defined as when the maximum density reaches 
10~^ g cm~^) for the radiation hydrodynamical calculations with 
/3 = 0.001 — 0.04. The lifetimes of the pre- stellar disc phase are all 
well converged for ^ 3 x 10^ particles, with the exception of the 
13 = 0.001 case. For the /3 = 0.005 - 0.04 cases, the lifetime of 
the first core/pre- stellar disc phase depends on the development of 
dynamical non-axisymmetric instabilities due to the high ratio of 
the rotational to gravitational potential energies. In particular, the 
transport of angular momentum via gravitational torques from the 
spiral density waves removes rotational support from the gas and 
drives the object to the second collapse phase much more quickly 
that it would evolve due to accretion and radiative cooling alone 
(see Section [T2| . However, the /3 = 0.001 case is rotationally sta- 
ble. Therefore, it evolves to the second collapse phase through ac- 
cretion and radiative cooling, but also through angular momentum 
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Figure 17. Snapshots of the column density viewed along the rotation axis at the point of stellar core formation (defined as when the maximum density reaches 
1 X 10~^ g cm~^. From top to bottom, the different rows are for molecular cloud cores with (3 = 0.001, 0.005, 0.01, 0.04. Note that the spatial scale is 
different for each row. From left-to-right, the columns are for calculations performed with 1 x 10^, 3 x 10^, 1 x 10^, and 3 x 10^ SPH particles. For each 
rotation rate, the evolution is qualitatively the same when using 3x10^ particles or more. 



transport by numerical shear viscosity. Numerical shear viscosity 
will transport angular momentum outwards, reducing the rotational 
support of the gas in the centre of the first core, driving it more 
quickly to higher densities and temperatures and, thus, hastening 
the onset of the second collapse. As the resolution is increased, the 
numerical shear viscosity decreases. Therefore, it is most likely that 
the longer lifetimes of the pre-stellar disc phase for the /3 = 0.001 
case with increasing resolution is due to the effect of the numeri- 
cal shear viscosity in SPH. We also note that the dominant term in 
the numerical viscosity (the term) is first-order, which explains 
why the convergence with increasing resolution is quite slow (left 
panels of Figure p^. Thus, for the first cores which do not undergo 



dynamical rotational instabilities, the lifetimes obtained in this pa- 
per should be treated as lower limits. 

In Fig. [it] we plot images of the first cores/pre-stellar discs at 
the time of stellar core formation for the /3 = 0.001 — 0.04 radi- 
ation hydrodynamical cases obtained using the different numerical 
resolutions. For the (3 — 0.001 case, the first core is right at the 
boundary of rotational stability, with a very weak spiral instability 
visible in the highest resolution case. For the /3 = 0.005 case, all 
resolutions except the lowest (1 x 10^ particles) resolve the rota- 
tional instability. As the resolution is increased, the sharpness of 
the spiral density waves improves and the pre-stellar disc becomes 
slightly smaller, but the structure definitely seems to be converging. 
For the /3 = 0.01 case, there are only two resolutions. Both pro- 




Figure 18. Snapshots of the column density viewed perpendicular to the rotation axis 25 years after stellar core formation (defined as when the maximum 
density reaches 1 x 10""^ g cm""^) in a /3 = 0.005 molecular cloud core. From left to right, the different panels are for radiation hydrodynamical calculations 
performed with 1 x 10^, 3 x 10^, 1 x 10^, and 3 x 10^ SPH particles. The bipolar outflow and the Shockwave propagating through the disc are similar in 
all cases, regardless of the resolution. 



duce a central object surrounded by a disc of radius ^ 90 AU which 
fragments. The higher resolution case forms two fragments, while 
the lower resolution case only has one fragment. Thus, while the 
number of fragments has not converged, the size of the pre- stellar 
disc and the fact that it does fragment seem to be robust. Finally, 
both resolutions of the /3 = 0.04 case produce ring-type first cores 
which fragment into four separate cores. However, with four ob- 
jects the subsequent evolution becomes chaotic. In the 3 x 10^ par- 
ticle calculation, the four fragments survive (at least until the first 
undergoes second collapse), but in the 1 x 10^ particle calculation 
two of the objects merge producing a triple system (see the bottom 
rows of Fig. [it] and |3]). 

Following stellar core formation, as discussed in detail in Sec- 
tion |3.3| a bipolar outflow is launched perpendicular to the disc and 
a Shockwave propagates outwards along the plane of the disc. These 
outflows and Shockwaves occur in all of the radiation hydrodynam- 
ical simulations, with the exception of the spherically-symmetric 
/3 = case. This result does not appear to depend on resolution 
(see Figure p^. However, as discussed in Section |3.3| and illus- 
trated in Fig.|14| the initial mass of the stellar core and its accretion 
rate are affected by resolution. 

In summary, we conclude that the lifetimes of the first 
cores/pre- stellar discs are well converged for initial conditions that 
produce rotationally-unstable pre-stellar discs, but that for more 
slowly-rotating cores (i.e. /3 = 0.0005 and /3 = 0.001) the life- 
times are lower-limits. We also find that the qualitative outcomes 
(rotationally- stable first core, pre-stellar disc with spiral arms, frag- 
menting pre-stellar disc, and torus) are numerically-converged but 
the details (e.g. how many fragments) may not be. 



